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A multilevel adaptive sparse grid stochastic collocation 
approach to the non-smooth forward propagation of 
uncertainty in discretized problems 

Robert L. Gates*^ and Maximilian R. Bittens*^ 


Abstract 

This work proposes a scheme for significantly reducing the computational complexity of discretized 
problems involving the non-smooth forward propagation of uncertainty by combining the adaptive 
hierarchical sparse grid stochastic collocation method (ALSGC) with a hierarchy of successively 
finer spatial discretizations (e.g. finite elements) of the underlying deterministic problem. To 
achieve this, we build strongly upon ideas from the Multilevel Monte Carlo method (MLMC), 
which represents a well-established technique for the reduction of computational complexity in 
problems affected by both deterministic and stochastic error contributions. The resulting approach 
is termed the Multilevel Adaptive Sparse Grid Collocation (MLASGC) method. Preliminary re¬ 
sults for a low-dimensional, non-smooth parametric ODE problem are promising: the proposed 
MLASGC method exhibits an error/cost-relation of e ~ ^-o.95 therefore significantly outper¬ 
forms the single-level ALSGC (s ~ and MLMC methods (e < 


1 Introduction 

It is well-known that many problems in engineering are subject to uncertainty in the input pa¬ 
rameters, e.g. material data, boundary conditions, and geometry. In the present case, we are 
interested in the forward-propagation of such uncertainty to the quantity of interest, e.g. deforma¬ 
tion or stress, which is usually obtained from the solution of a partial differential equation. Three 
classical categories of methods exist for the solution of such stochastic partial differential equa¬ 
tions, namely the Stochastic Galerkin (SC) [1], Stochastic Collocation (SC) [2], and Monte Carlo 
(MC) [3] approaches. The selection of an appropriate method depends strongly on the number 
and independence assumptions of the random variables, as well as on the smoothness of the solu¬ 
tion in stochastic state space. Regardless of the employed method, stochastic partial differential 
equations pose a challenge to being solved efficiently when the underlying deterministic problem 
is computationally costly. It is the purpose of this work to contribute to this field by proposing 
a method for reducing the computational complexity for a wide range of moderate-dimensional 
parametric problems encountered in engineering practice. 

In the following, we consider a general class of problems from the field of computational mechan¬ 
ics, where solutions to inequality constrained stochastic partial differential equations are sought. 
Prominent examples include the computation of complex material behavior with uncertain ma- 
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terial parameters as well as the contact of deformable bodies with rough surfaces in a generally 
nonlinear setting. Due to the complimentary condition, such problems are classified as being 
non-smooth in the random parameter domain [4]. In both exemplary cases, the rough surface as 
well as the uncertain material parameter field is represented by a stochastic process which, pro¬ 
vided that it has bounded second moment, can be approximated by a truncated Karhunen-Loeve 
expansion. This leads to the assumption of a multilinear combination of independent random 
variables parametrizing the deterministic problem. Due to independence, it is possible to utilize 
a double-orthogonal polynomial basis of tensor-product structure for the state space, leading to a 
decoupling of the random dimensions (see e.g. [2,3,5]). This choice of basis renders the SG and 
SC approaches equivalent, resulting in a non-intrusive method. 

Unlike Monte Carlo integration, the convergence of classical SG and SC methods relies on regu¬ 
larity properties of the quantity of interest in the random parameter domain [2,6]. Unfortunately, 
the class of problems considered herein clearly violates the required smoothness assumptions [4], 
making MC methods an attractive alternative, despite their slow convergence. However, in many 
cases, areas of reduced regularity are confined to certain regions of state space, suggesting an 
error-adaptive approach to sparse grid SC methods [7-9] for problems involving a moderately 
large number of stochastic dimensions. Such adaptive methods utilize local hierarchical basis 
functions with tensor product structure, naturally providing for an improvement estimate for 
each adaptively refined collocation point while overcoming the oscillations incurred by the use of 
global interpolating polynomials. 

This work proposes a scheme for significantly reducing the computational complexity of dis¬ 
cretized problems involving the non-smooth forward propagation of uncertainty by combining 
the adaptive hierarchical sparse grid stochastic collocation method [7-9] with a hierarchy of suc¬ 
cessively finer spatial discretizations (e.g. finite elements) of the underlying deterministic prob¬ 
lem. To achieve this, we build strongly upon ideas from the Multilevel Monte Carlo method 
(MLMC) [10-13], which represents a well-established technique for the reduction of computa¬ 
tional complexity in problems affected by both deterministic and stochastic error contributions. 
The resulting approach is termed the Multilevel Adaptive Sparse Grid Collocation (MLASGC) 
method. It is remarked that previous works on the topic of multilevel methods in adaptive sparse 
grid stochastic collocation methods [14] focus on the acceleration of iterative solvers by using a 
low-fidelity interpolant of the stochastic state space as an initial guess for newly added collocation 
points. It is emphasized that, while the underlying ideas are very similar, our approach is more 
classical in that we consider the term “multilevel” to apply to a hierarchy of successively finer 
spatial discretizations, as suggested in [15]. 


2 The multilevel adaptive sparse grid collocation method 


We begin by summarizing the so-called Adaptive Lagrangian Sparse Grid Collocation method 
(ALSGC) [7,8,16] and subsequently extend it to the use of multilevel deterministic discretizations. 


2.1 Construction of an initial sparse grid 

For the construction of a d-dimensional sparse grid on U x ... x x ... x by the Smolyak 
algorithm, introduce the multi-index i = {ii,... ,ij,... ,id), in which ij € N+. Further, define 
Si = {i I 1 — d-|-|i| = 1} as the set of multi-indexes belonging to level I of the sparse grid. In order 
to relate the level ij of the j-th dimension to a certain number of univariate points, we choose a 
nested rule, e.g. the Clenshaw-Curtis rule: 


n{ij) 


1 if ij = 1 

2*-i-gl ifi^>l 


( 1 ) 
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and denote 0^^. as the resulting set of n{ij) univariate points on the interval Ij. For each level I, 
the sparse grid is constructed by the relation 

= U X ■ • ■ X X ... X ©iJ . ( 2 ) 

The adaptive grid discussed in the following section requires an initial set of sparse grid points 
up to level iinit as a starting point for refinement. Hence, we define the initial grid ©init as the 
union of the incremental, i.e. level-wise, grids 0/: 

L[nit 

©init = U ■ (3) 

Z=1 


2.2 Tree structure of sparse grid points 


Due to the nested structure of the collocation points, the sparse grid introduced in the previous 
section admits a k-aiy tree structure. In particular, every parent collocation point on level I 
has at most two children per dimension on level I + 1, i.e. d < k < 2d. In order to make 
this notion precise, we identify univariate points on level ij of the j-th dimension by integers 
rrij’ G {m G N’*' : m < n{ij)}. Given a univariate parent collocation point on level ij identified 
by the index , its index on level ij -I- 1 is given by 




2 if ij = 1 

2nij — 1 liij > 1. 


(4) 


It is emphasized that an existing collocation point identified by the index nij exists only on 

level ij, while its index nij on level ij + 1 simply corresponds to a non-existent placeholder 
resulting from the nested structure. This convention allows for the identification of the at most 
two univariate children on level ij -|- 1 to point rrij as 

= I - l} if + f) (5) 

y ^nij~^^ — 1, 4" l} if otherwise . 

For a univariate point identified by the integer nij as well as the level ij, it is straightforward 
to recover its coordinate. For instance, in the case of univariate points distributed evenly on the 
interval Ij = [—1,1], the coordinate is 



-1 




if n{ij) = 1 
if otherwise. 


( 6 ) 


Turning back to the multi-dimensional case, collocation points on level I are identified by the multi¬ 
index m* = ,..., nij , ■ ■ ■, rn^), where i G Si. We are then able to identify its coordinates 

^(m*) = X ••• X X ••• X ||(to^‘')| (7) 

as well as the indices of its children 

cr = {{mi} X ... X {rn/Sl] x x x ... x {m;/} | j = l,...,d|. (8) 
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Finally, we denote 


Cl 


{{{(in}} if; = i 

{ct" I m* eCi_i} if/>l 


(9) 


as the set of all collocation point indices on level 1. For an example of the hierarchical construction 
in two dimensions, see tab. 1. 


2.3 Local hierarchical Lagrange interpolation 


As already mentioned, the ALSGC approach uses local hierarchical basis functions to overcome 
the drawback of oscillations, well-known from interpolations using a global Lagrange polynomial 
basis. We begin with the construction of univariate basis functions, which will then be extended 
to the multi-dimensional case. 

For identifying the support of the hierarchical basis, every point index rnj of an existing 
collocation point on level ij > 1 possesses at most two neighbor point indices 

|m}'’ -1-11 ifm}"=l 

|m*'’ ~ l| ifm}^ = n(ij) (10) 

— 1, rrij + l| if otherwise . 

per dimension j. It is remarked that these neighbor point indices need not correspond to existing 
points. In case of non-existence, they correspond to placeholders in the nested structure. We 
emphasize that the linear basis function with significant points 



Pi/ = 


l^'(n) I n e n™" I U ) } 


( 11 ) 


has bounded support supp^a™^^ 


• r pyii Pi 

mfp-d, suppj 

C(e) = 


and fulfills the following conditions: 


1 if^ = ^(m*J 

0 if ^ ^ supp(o™), 


( 12 ) 


The choice of linear basis functions is made for reasons of simplicity. Other possible candidates 
include multi-resolution wavelet basis functions [9] and higher-order Lagrange polynomials [17]. 
In all cases, the basis function for level ij = 1 is defined as 

a{ = 1. (13) 


Returning to the multidimensional case, the d-dimensional basis function can be constructed using 
tensor products 


d 


< = 


O a, d O ■ 




(14) 


i=i 

of univariate basis functions. The set of all d-dimensional basis functions at” on level I is denoted 
by Ai = {at” I m® € C;}. For an example of the hierarchical basis in two dimensions, see table 2. 
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The use of hierarchical basis functions subdivides the sparse grid interpolation space 

L 

Vt = ®Wu (15) 

1=1 

into an orthogonal sum of hierarchical difference spaces 

Wi = spanja^^ : € Ai} . (16) 

The interpolant of the hierarchical difference space for a function / : > R is defined as 


f E «r( 0 -/(^M) 


if; = 1 


E «r( 0 - 

\ Cl 


-Zz_i(/)(^(m®)) 

if; > 1 

(17) 

= Y. aT{i)-wT. 



(18) 


Cl 


where ic™ denotes the hierarchical surplus belonging to the basis function on level 1. The 
hierarchical surplus on level / = 1 is simply the function value at the coordinates of the collo¬ 
cation point located at = {ml,... ,mj,... ,m;)) = (1,..., 1,..., 1)). For levels I > 1, the 

hierarchical surplus is constructed using the difference of the function value at the coordinates of 
a collocation point ^(m®) where i G Si and the value of the interpolation Z;_i of level Z — 1 at 
the same coordinates. The complete interpolant is then recovered by the sum over all hierarchical 
difference interpolants: 


L 


1=1 


(19) 


2.4 Adaptive refinement 

The hierarchical surplus provides for an improvement prediction of the local interpolation when 
transitioning from level / to / -|- 1. It therefore comes naturally to use the hierarchical surplus as 
the criterion for the adaptive refinement. Hence, for a given tolerance e^ef € R"*", children c™ of 
a collocation point on level I identified by the multi-index m® are created if \wY^\ > e^et- This 
allows for the identification of the set of adaptively refined collocation points on level / -f 1 as 

Cz+i = {c7®|m®eC,A|u;ri>eref}. (20) 

It is remarked, that for a proper adaptive refinement the initial level Tinit of the sparse grid 
©init must be chosen sufficiently high. Otherwise, the possibility increases that the refinement 
stops prematurely when a function value at a refined point coincidentally equals the value of the 
previous level interpolant. 


2.5 Multilevel splitting 

In the most general case, it is assumed that the quantity of interest u{-,^) of the underlying 
deterministic problem can be approximated to a certain precision level r by an arbitrary scheme, 
e.g. an ordinary differential equation integrated using a difference scheme with time step Atr, a 
partial differential equation solved via finite elements with characteristic mesh-width hr, or simply 
a 16 • 2’’-bit floating point precision evaluation of a function. In practice, we may be interested in 
computing the underlying deterministic problem to a precision level R, hence it is straightforward 
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Figure 1: Exact solution u{t = 1,^) of ODE G [—1,1]^) 


to acknowledge that the telescopic sum 


R 

U_R = Ml + Ur — Ur-1 (21) 

r=2 

realizes this requirement. It is remarked that computations to precision r are usually significantly 
more costly than computations to lesser precision r — 1. 

Due to the hierarchy of successively finer discretizations, we assume that the variance of the 
level correction V[ur — Mr-i] —> 0 as r —>■ oo. In analogy to MLMC methods, it appears sensible 
that the number collocation points required for achieving a given interpolation error tolerance 
is related to the variance of the interpolated function. Hence, the variance decay of the level 
correction Ur — Ur-i is exploited in order to reduce the overall cost of the computation. In 
particular, we estimate ui as well as the subsequent corrections Ur — Ur-i independently using 
the ALSGC method. We then recover the response surface approximating ur{^) by 

R 

I^^[urU) := I^^luim + Y.I^^[Ur - Ur-iKO ■ (22) 

r^2 

Subsequent integration of the response surface X^^[t6i^](^) for purposes of stochastic moment 
estimation then follows in a straightforward manner. 
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3 Numerical results 


For a preliminary numerical investigation, we choose a parametric € [—1,1]^) first-order linear 
ordinary differential equation 

+ (|2 - (Cl - 1)^ - (C 2 - 1)^1 + S)u = 1 (23) 

with initial condition u{t = 0,C) = 0 as well as a regularization parameter S = 10“^. For reasons 
of simplicity, we will only be interested in the final value u(t = 1,C)- The problem admits an 
exact solution 


u{t,0 


1 - exp(-t(|2 - (Cl - 1)2 - (C 2 - 1)^1 + S)) 
|2-(Ci-i)2-(C2-i)2| + ^ 


(24) 


which is shown in fig. 1 and used as a reference for the numerical solution Ur(t = 1,C)) obtained 
using a forward Euler integration scheme with time-step Atr = 30“^ • 2“’’. 

We define the total error, including interpolation as well as discretization contributions, in the 
T^-norm as 




L2([_i,i]2) 



w(C)-TML/sl[^^](^) 



1 

2 


(25) 


Figure 2 shows this error, computed by Monte Carlo integration using 10® points, achieved by 
the MLASGC and ALSGC methods over computation time. The data points correspond to 
discretization levels i? = 4,..., 15. The adaptive refinement tolerance Cref = AIr for the single- 
level ALSGC method is chosen equal to the global truncation error Cdiscr ~ Atu of the forward 
Euler method, aiming to balance the contributions of discretization and interpolation error. For 
the MLASGC method we distribute the desired overall refinement tolerance Cref = Atn across the 
R multilevel interpolants in a linear manner such that 


^ref 


R 



r=l 


2reref 


(26) 


The linear increase of the refinement tolerance across levels r aims to enforce a high refinement 
tolerance on coarse discretizations and a lesser refinement tolerance on fine discretizations where 
additional collocation points would be costly. This choice performed better than the uniform 
tolerance distribution = e^ef/R- The idea behind this rather heuristic attempt to balancing 
cost is confirmed by the number of points required per level r of the multilevel interpolant for 
achieving an overall refinement error in the order of Atn (see fig. 3). It is however emphasized 
that even a uniform distribution of the refinement tolerance leads to a significant decay of the 
required number of collocation points across level corrections, such that the linear distribution 
should only be considered a slight correction. 


4 Conclusions 


The preliminary results for the low-dimensional, non-smooth parametric ODE problem considered 
herein are promising: the proposed MLASGG method exhibits an error/cost-relation of e ~ ^-0-95 
and therefore significantly outperforms the single-level ALSGG (e ~ and MLMG methods 

(e < ® [4,18]). Due to a lack of mathematical analysis of the new MLASGG method, no 

special cost/error-balancing is performed, leaving room for further optimization. It remains to be 
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investigated if the techniques presented in [15,19] in terms of non-adaptive multilevel collocation 
methods are similarly applicable to the MLASGC method. It is also emphasized that the new 
method is not limited to the use of the ALSGC method for the interpolation of the level correction. 
In fact, a further performance increase would be obtainable by the use of a multi-resolution 
hierarchical wavelet basis in state space which introduces a true local error estimate due to the 
fulfillment of the Riesz property (see e.g. [9]). Finally, we remark that the non-intrusive nature that 
the MLASGG method shares with other variants of collocation methods lends itself excellently to 
parallelization and implementation into existing deterministic frameworks. 
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5 Appendix 

5.1 Example: Grid and basis function construction in 2D 


1 

index-sets Si and children-sets Ci 


1 

Si = {i\ii + 12 = 2} 

= {(1,1)} 

G = {{{(1,1)}}} 



L 

'(1.1)(13) 

2 

^2 = {i\ii+i2 = 3} 

= {(2,1),(1,2)} 

= {{1,3} X {!},{!} X {1,3}}} 

= {{{(1,1), (3, l)}-(2d), {(1,1), (1,3)}*=(b2)}} 

< 

• (1.3)(1>2) 

< 

(3.1)P.i)' 

3 

^3 = {i\ii+i2 = 4} 

= {(3,1),(2,2),(1,3)} 

n _ L(l.l) (3.1) (1.1) (PS)! 

^3 — 1 ‘'(2.1)>‘'(2.1),‘'(1.2),‘'(1,2)J 

= ({{2}x{l},{l}x{l,3}}, 

{{4}x{l},{3}x{l,3}}, 
{{l,3}x{l},{l}x{2}}, 
{{l,3}x{3},{l}x{4}}} 

= {{{(2,1)}(3.i),{(1,1),(1,3)}(2.2)}, 

{{(4,1)}(34),{(3,1),(3,3)}(2'2)}^ 

{{(1,1),(3,1)}(2.2),{(i,2)}(i.3)}, 

{{(1,3),(3,3)}(2.2),{(i,4)}(i.3)}} 

• 

(1.3)(2.2) 

(1.4)(i.3), 

(3.3)P.2)* 

1 

(4,1)(3.i) 

k A A 

< 

(1,1)(2,2) 

• 

J W W 

• (1.2)(1’3) 

(3.1)P.2) 

• 


Table 1: Example: Construction of the sparse grid level 0;, index-sets Si and children-sets C; 
2D 
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